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ABSTRACT 


The technique of reconstruction tomography used 
in radiology is applied to measure density in bubbly air- 
water flows. The ‘Convolution 1 algorithm used, converts the 
fan beam data to a parallel beam data and then a suitable 
filter is used. (in our case the generalise filter). Then 
convolution is carried out, followed by superimposition. 
Printing of density values is done in the format of density 
maps. The Generalised filter used proved superior to the 
Shepp and Logan filter which was used as a reference. 

The error band in the measured values of density 
is -0.03 and its underestimation is fairly constant and 
predictable within the range of 0.6 to 1.0 . 

The New Algorithm is an alternate method of 
r e.'.constructing densities but the C.P.U time is more. The 
error was practically zero, for a certain class of density 
distribution called radially symmetric distribution. 
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CHAPTER 1 

INTRODUCTION 

One of the major probl ens confronting a Nuclear 
reactor safety engineer, , is the measurement of void fraction 
or density during Loss-Odf-Coolant-Accident (LOCA). The 
development of an accurate and reliable two-phase flow 
measurement technique for various flow systems and components in 
light water reactors is of great importance and significance. 

The present study looks at such a new vista, which will 
alleviate the problems in the chemical industry, food processing, 
waste treatment and fludized flow systems, including fludized 
bed combustors. The method used is that of computerised axial 
tomography. 

The most familiar images are those formed directly 
by optical instruments using visible light reflected or 
transmitted by an object. However, in many applications 
in which an image is required, we can make only indirect 
measurements by probing the object with invisible radiation 
or by interpreting it. Often the measured data is not in a 
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form suitable for immediate interpretation, but is related 
to the image in a known way. The general aim of all image 
reconstruction procedures or algorithms is to process the 
data to form an image and So to facilitate the interpretation 
of measurements. 

To determine the density of the object under test 
various strip integrals corresponding to a particular angle of 
view were taken and the set of strip integrals are' called 
projection of the object. Given a nomber of such projections 
at different angles of view, the estimation of the corresponding 
distribution within the object is the basic problem of image 
reconstruction from projections. Computerized Axial Tomography 
(CAT) is the most significant application to date of image 
reconstruction from projections. 

According to Webster, 'TOMOS' means cutting and 
graphy is study and together mean a special X-ray/Gamma ray 
technique by which detailed images of a structure lying in 
a specific layer of tissue may be obtained, with images of 
structures in other layers eliminated or blurred. 
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A basic and crude form of CT (computerized 
tomography) was practised in Japan as early as 1947 
and called as Rotatography by taking several X«*ray images 
at different angles { .0-360°) around the patients body. 

The practice now days is to use a detector to 
collect the data and process using a computer/ which has 
the capability to detect minute differences in X-ray/Gamma ray 
attenuation coefficients, instead of directly using X-ray 
pictures (Radiograms) as was done with Rotatography . 

The present study involves "Convolution" algorithm. 
This algorithm is considered as highly efficient by numerous 
investigators in computerised axial tomography. The 
reconstruction method in use in all CTs is the Convolution 
method . 

The above method has been aptly been demonstrated 
by Kulacki et al (1). Reconstruction tomography was 
adapted to the needs of two-phase flow density measurements. 
Design of a prototype scanning densitometer has also been 
suggested by Schlosser et al (2). 
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There are several parameters in the Convolution 
algorithms and becomes implementation dependent. The present 
study considered all the free parameters and their effect on 
the problem in question. The parameters considered are viz:- 

1 . Number of scans =- N - 

2 . Number of beams in a scan - M - and 

3. Various filter functions 

The reconstruction algorithm used is that of Kwoh 
et al (3) with the Shepp and Logan filter. The present study 
investigates other possible filters like Ramachandran filter. 
Parabolic filter and Generalised filter. 

The present study has been undertaken with the main 
objective to improve the reconstruction algorithm used in the 
work of Kulacki et al (1). Especially the filter function 
( Shepp and Logan ) used has been carefully investigated. 

A FORTRAN program for the algorithm has been written 
and implemented with flexibility to change number of Rays per 
view and number of projections in scan as per the computer 
time and accuracy needed. 



Specific objectives of the study are:- 


1* To find the optimum values for number of Rays 

per projection (LITM). 

2. To find the optimum values for number of projections 
in scan (St) and 

3. To find the best filter.. 

LITM was checked for two values 20 and 40 . 

N was checked for three values 10, 20 and 40. 
Parabolic filter, Ramachandran filter and four cases 
of Generalized filter were investigated, graphs plotted. In 
this test for best possible filter Calibration curves were 
drawn for a few known densities of material like Air, Water, 
Walnut and Pine. 

The data used in all the computations have been 
taken from Kulacki et, al (1), 

A 'New Algorithm' was investigated to get exact 

reconstruction. This was developed by Prof .R .K .S.Rathore. 

The data used were test values. Experimental Data is to be 
adopted for this metjjod. 
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CHAPTER-- 2 

LITERATURE REVIEW 

2 . 1 OVER VIEW 

One of the complex situations to analyse in 
Nuclear reactor safety is the two-phase flow condition 
and void fractions formed in them* There are various methods 
used for this purpose and Kulacki et al (1) have amply and 
sufficiently proved in their study the use of scanning 
gamma-ray densitometer. There are no disturbance caused 
which are inevitable if probes are inserted in the flow for 
point density measurements. Also a scanning gamma-ray 
densitometer offers the potential for obtaining both point 
and cross-section averaged values of density without disturbing 
the flow field ( because of invasive technique ). Most of the 
present study is based on the technique involved and developed 
by ( 1 ) . 

The other prominent in the category of invasive 
type of measurement of void fraction are :~ 

Nuclear Magnetic Resonance ( NMR ) and 


Ulltrasonics ; 
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Computerised axial tomography (CAT) the method 
adapted to two— phase flow in Nuclear reactors by Schlosser 
et al (2) with the aid of Kwoh et al (3) algorithm is 
explained in the following pages. 

Ramachandran et al (5) has explained in their 
study a simplification of reconstruction by using 
convolution instead of Fourier transform, which greatly 
reduces the computer time and space needed alongst with a 
better accuracy. 

2 . 2 RAMACHANDRAN FILTER 

This new technique is used for 3D reconstruction 
from their transmitted radiographs. Convolution is used 
in the real space of the object, without using Fourier 
Transforms (F.Tts). The object is rotated about an axis 
at right angles to the 'direction of a parallel beam of 
radiation, and sections of it normal to the axis are 
reconstructed from Data obtained by scanning the corresponding 
linear strips in the shadowgraphs at different angular 


sett ings 
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Since convolution method involves only one 
summation over one variable at a time, while F.T requires 
two, convolution is much faster. It is also shown by tests 
that convolution method is more accurate* 


The basic step in the solution of reconstruction 
problem consists of replacing the three dimensional 
density distribution by a set of two-dimensional density 
functions in a series of sections perpendicular to . 
Z-axis. 


oo 

r 

F(Rj.e) = ' g(l;0) exp( 2 w iRl ) dl (1) 

f(r;0) = J ! F(R,0) exp (-2^'iRr Cos(0-£))R dR d© 
0 f) ....(2) 

Suppose we define: 

oO 

g*( 1 , 0 ) = j IrJ F(R;8) exp (-27T iRl)dR (3) 

JsO 

Equation (2) for f(r,<$) becomes 

T 

fir, ft) = f g (r Cos {$-&),&) (4) 

o 

Fourier inverting (1) 

r 

g(lj&) = ! F(R;0 ) exp(-2 7f iRl )dR (5) 

-GO 

F.T of g(l;0 ) = (F.T. of g(l;9)) x (F.T. of q (1)) 

R is the F.T. of the function q (1). 

, f Rj = / q(l ) exp(27TiRl) dl (6) 

-cO 
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Using convolution theorem for the inverse of the product 
of F.T's 

g 1 (1; 0 ) = y g(l 1 ;^) .qd-^) dl 1 ......(7) 

^’=Oq 

q (1) = - f /R/ exp(-27T iRl ) dR ......(8) 

_y»o 

The integral diverges for the limits so replacing the 
limits by -A/2 and +A/2 where A= 1/a. 

AL. 

q A (l)= / (R j exp (—2 Tf iRl ) dR (9) 

q(na)= l/4a for n=0 

= -1 /fT~ h 2 * a^" for n odd (10) 

= O for n even 

f(r ,Q) is the two dimensional density distribution. 

fCr^^rs) will be the data of three dimensional 
density distribution using all linear strips at Z in the 
different shadowgraphs taking at various angles O . 

F(R;©) is the Fourier transform of g(l;9) 
where g(lj$jz) is the logarithm of the ratio of intensity 
at point P(z,l) to the incident intensity . (which is a 
measure of the linear integral). 

The biggest problem in this filter is that it is 
zero for all even values so the filter becomes highly 
problem specific. 
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2.3 SHEPP AND LOGAN FILTER 

The study of shepp et al (6) .reveals that the 
Fourier reconstruction may be viewed simply in the spatial 
domain as the sum of each line integral times a weighting 
function of the distance from the line to the point of 
reconstruction. 

Ramachandran’s linear interpolation is merely 
the choice of a particular weighting function. A modified 
weighting function simultaneously achieves accuracy, 
simplicity, low computation time as well as low sensitivity 
to noise. 

Using a simulated phantom, comparison studies were 
done on the Fourier algorithm and a search algorithm very 
similar but notidentical to Housfields algorithm. The search 
algorithm requires 12 iterations to obtain a reconstruction 
of accuracy and resolution comparable to that of Fourier „ 
reconstruction and more sensitive to noise. To speed the 
search algorithm using fewer iterations results in error. 

The steps similar to Ramachandran 1 s filter 
derivation are done with the modification in the weighting 


function 
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q(0) = 4/77 a 2 " for K=0 

q (Ka >=~4 /ff a~ ( 4K^-1 ) for K=+ 1,2, 

This is the Shepp and Logan filter, 

2.4 GENERALISED FILTER 

Jn their study Kwoh et al (4) have given a new 
filter called generalised filter. One general class of 
digital w - filter without interpolation is defined by 


W (w ) = 


/ 1 j I _| W ^ 1 - A j wf ^ ^ IT/ 

■ — w exp — ? = — w exp w 1 w 

a 1 j a f ' 


1 ! 2TT 
a J a 


»w exp , 

t J ( a 


ilp; . J p 


<« <; 2jr 
a 


where zeta ( <^ ), p^ O? a is the sample spacing and W(w) 
J 

is the periodic with period 2 TT/a . 

il P 


Zeta = «, 


I W d 


is called the damping factor. It is a 


function of the cut-off frequency w c and the parameter p. 
the increased number of parameters over what was used in 
other filters makes it easier to fit a filter specification. 

tw >P 

The portion exp ( - j — j ) serves as a -"switch. 

' ct 

At low frequencies , it is not triggered and W(w)= jw J . 

For higher frequencies, it switches to zero as w approaches 
the cut-off frequency w c . W(w) tends to zero with a rate 
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depending on the parameters p and w^. To minimise oscillations 

it is desirable to choose w and p so that a smooth cut off 

c 

results. This is exactly what that is done. Four such values 
taken and tested to find which amongst them fit in and give 
the best of results for the problem in question. 

When ^ =0, the generalised filter is identical 
to the triangular wave filter used by Ramachandran. For certain 
values of ^ and p, the generalised filter approximated 
the filter of sheppand logan very (fairly) accurately. 

For p>l as in all the four cases considered numerical 
analysis is used to obtain the response. 


W(na ) = — 


/A 

wexp- 


w 1 


Cos naw dw 


( 11 ) 


This is the response for n = O 


for other n 1 s , We change variables by putting z=naw. 


So that (11) approximates to 


w 


(na)= 




rtk^j z ^ 


J z p 

_c(“ n ^ Cos 2 for n=l/2,.. 


( 12 ) 


for n=l we get, 

1 C ^ • } -7 P 

W(a) / z exp_J(~ ) Cos z dz 


( 13 ) 
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for n>l, (12) is put in alternate form 


when n is even 

i -4(z+27Tj) 

w(na)= jj-g-'rf-j Cos z ( ;> (z+27f j) exp9(<2rf) 


1) p 

dz ) 


w 


j-c 


b 

hen n is odd > 1 
W( na)=„ a i^j * ( ^lz+2 ZT j) ®P d 2 

' — L> 

-^< z+ 2y(n-l) P az) 


,7T 


Cos z ( (Z+2 7T(^“ 1 )) exp 

•s/ 


0 n 


W(na) = W(-na) for n= 1,2 

This filter in turn was used with various values 
for its parameters p and Zeta, •‘.which were recommended by 
Kwoh et al . (4 ) . 

P= 2.03399 and Zeta = 2.243E=06 is the approximation 
of a generalised filter to Shepp and Logan filter. 

P= 1.2 and Zeta=5E-04 is the approximation between 
Shepp and Logan and parabolic filter. 

P=2.0 and Zeta=3E-05 best minimized the internal 
undershoot while preserving the shape of the object. 


P=2.6 and Zeta=0.96E-05 was also used 
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2.5 PARABOLIC FILTER 


The parabolic filter is an approximation of 


certain parameters of a generalised filter and particular 
value of p and Zeta ( ). 


The filter has response of 

7T 


W(0)= 

W(na)= 


3a ^ 

-1 

Jf ( na)‘ 


for n=0 


for n= + 1 7 +2 


w 2 / 


W(W>= i| H- 06 wi f 


is the frequency response of a parabolic filter. 

2.6 Filter used : 

Three dimensional (3D) reconstruction for diverging 
X-ray/Gamma ray beams as studied by Kwoh et al (7)has shown 
that fast digital convolution technique is extended to 
diverging X-ray/Gamma ray beams by suitable changes to parallel 
X-ray/Gamma ray beams conditions. This along with the 
generalised (wj filter yields good accuracy and excellent', 
results* It also has a flexibility to cope with noise and 
substantial reduction in X-ray dosage. 


. t 
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The background far this work is based on 
Bracewell and Riddle who first suggested Fourier Transform 
( F.T ) might be used to reconstruct an object and that X-ray 
scanning replaced 3D object by a series of 2D cross-section 
of the object. Then Ramachandran and Laxminarayan * s work , 
and their significance already discussed have improved the 
situation further by introducing convolution. 

The key equations are 


f(x,y) =/^ Tr v4J w i dw exp(iw(xCos^+y^Sj.T)f£) ) .(14) 


7 2 iry 


•‘-OO 


w 


is spatial frequency and p(w,e?) is F.T. of 


projection data for parallel beams. 

b(t,&) = exp( (wt ) ) dw 


..(15) 


t= xCos{H-ySin 8 , b(t,£)) is called the back projection 


formula. By the convolution property of the F.T 


B) w(t-^)d^ ...(16) 

where p(t / 6') and w(t) are respectively inverse 'F.T of 
p(WzB') and/w/. 

The inverse F.T of jwf does not converge so (16) 
is valid only in a limiting sense. Since most of the objects 
have spatial frequency of interest only belcw some cut-off 


b(t,0O= p( 



frequency we can replace w with a function w fw) 
which approximate |w| for w<w and zero for w^-w . The 
class of functions W(w) with these properties are called 
^w|functions. This is exactly the generalised filter functi 
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CHAPTER— 3 

PROGRAM DESCRIPTION 

3.1 Overview: 

The application of the concept of computerized 
tomography scanning in two-phase flow problems results 
in absolute point dens ity values in "the pipe cross section 
averaged over time. The program involves four basic steps 
which include data input, adjustment of data, reduction 
and finally the reconstruction of the density map. 

The data were taken from Kulacki et al (1) and 
was in the form of a fan beam geometry. Data was taken at 
every 2.5 degrees for a total of 50 degrees. The normalised 
counts were denoted by h( ) , refer Fig. 3.1 

where, A 

p= ioi for O ^ iv<C N 

\ _ — M 

A= Sa for ^ , 

M is the number of diverging beams in the fan 

N is the number of scans 

a ~ ^max /M, 
and o£= 1T/N 
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Since the tomography concept is applicable only to 
the parallel b6am case/ the fan beam data has to be 
transformed. The transformed input/ p(L,9 ) was determined 


by 


X = 1 0 = 9 -Sin“ 1 (L/D)/ 3.1 

r 4 " ' o { 

1— (L/D ) 


V 


and p(l/0 ) =h (\fi)= h (- ... ...... 

J 1— (L/D ) 


1 - . .l / © -Sinl ^ 

2 


Sint g)) 3.2 


equation (3.2) involves the lengths L,D and A , and angles 
© and/3, which are define! in Fig. (3.1). 

The third step involves the use of a filter 
function/ W and the computation of the convolution C (la). 


C (La ) = 


W 2 

<■ 

K= 


-M 

2 


p(t R ^ ) W( t L“ t k) 


where/ 

a ^h 

^i is the j angular position 

t^ = ia for i= j, k,L, 

j,k/L are all integers/ 

t .= X cos 6 i + Y sin^. 

3 r- 

Once C( t, Q . ) is known, the back projection, b and 

g x 


can be determined. 
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b (t,£K) = a 


+* -»t 

^L+l L 


C(t T A )+ t “ t L 


b/i + t L t C(t L+ 1 A) 

t L+l" t L L+1 1 


and finally the reconstructed density function f(x, y) is 


f(x,y) = Is ^ 


N-l 


i=0 


b(t, 9^}, 0$i^N-l 


The print out of the program is enclosed in 


Appendix A 


3.2 Steps involved 


The essential steps involved were as follows J 
The input data was given in the proper frsrmat along with 
source to center distance (D), diameter of the object (T) 
Number of scans and number of projections in scan (M), etc. 
Then this was supplemented with LITM and N and a conversion 
to parallel beam done, along with interpolation to increase 
the number of data points. 

Then the next step was in choosing the filter 
function. Schlosser et al (2) in their study have amply 
proved that Shepp and Logan filter is simple and good 
enough to the data used. This is verified and Parbolic 
filter and four cases of generalised filter are also taken up. 
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The generalised filter was calculated using 
integration routines and this was done with a limited 
accuracy. There were four integrations to perform to 
calculate separately values for N=0 , 1 , add and even. 

This shows as to why generalised filter was not 
used by Kulacki et al (1). The calculation of generalised 
filter is really a cumbersome process involving four 
integrations. Four types of generalised filters were 
considered depending upon which other common types ' 
then approximated to. There were four values of p and ZETA 
chosen and they are as follows. 

1. P= 2.0 and Zeta= 1.07E-02 

2. P= 2.0 and Zeta= 2££'0 2, 

3. P= 2.0 and Zeta= 4.75E-03 
4 . P— 2.0 and Zeta = 2 . 67E , “*03 

These values of filter were used and Convolution 
found. Then Superimposition and Printing done in the proper 
format and LITF directly gives density values. 



3.4 Data used : 


The Data was taken from the study conducted by 
Kulacki et al(l). The source of nuclear radiation was 13.86 mCi 
of Cs-137 and fixed. The Collimator fixed at the other end was 
rotated through and angle of 50° (25° on each side). 

The Detector was connected to this and the 
readings taken through preamplifier, amplifier, baseline 
restorer. Single channel analyser and a timer. 

The Data was collected for eight various materials. 
Four of them water, air, walnut and pine are used for 
calibration .purpose. The other four 90% density, 80% density, 
70% density and 60% density were taken for six types of 
scans and calculated and drawn accordingly on the graph. 
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CHAPTER-4 

RESULTS 

4.1 THEORITICAL CONSIDERATIONS 

Density measurements in two-phase flow via- the 
techniques used in axial computerized Tomography involves 
the data output format, interpretation of output and 
calibration of the computer algorithm. In each area the 
problem of relating what the reconstruction calculation means 
both physically and mathematically must be addressed. The 
data obtained in the present study is interpreted and best 
utilized for the overall objective of the investigation. 

INTERPRETATION OF DATA OUTPUT 

A typical reconstruction output appears in 
Matrix form shown in Fig. 4 . a The elements of this 
output Matrix vary from 0 to 99. At the bottom of the 
picture are indicated the maximum and minimum values of -•LITF* . 
This variable, LITF, is related to point density values. 

In fact, each number corresponds to the density of unit 
cell around it. Linear interpolution- is used to reduce 
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the actual values of LITF within a range of 0 to 99. 

Since the quantify of interest is the horizontally 
averaged density/ fy, the cross-sectional average of 
LITF is required. Thus ^LITF^- is the sum of the values of LITF 
over the cross-section divided by the number of unit cells 
present. 


The algorithm output is a 4Gx40 Matrix while the/ 
pipe wall effects prevents any reasonable values to be had 


above 31 units. Hence 31 x 31 Matrix is considered good and 
used for averaging and finding ’LITF 1 . 


%\ 

LITF =£ LITF 
*=*961 


.. . (4.1) 


where 961 is (31x31) the total number of unit cells. 

The reconstruction accuracy for Shepp and Logan 
filter is around + 10% . The accuracy for generalized filter 
is around + 4%. 


CALIBRATION 


The quantity ^LITF/determined by equation (4.1) 
corresponds to the actual density values across the flcwarea. 
Scans for the cases of Water Qp=1.0 g/cm ), Air (gp^approximately 
zero g/cm ) # Walnut (p=0.73 g/cm ) and Pine (jp=0.41 g/cm ) 
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result in four different values of LITF for known 
densities- Since these known densities are accurate within 
+ 0.5% a reasonably good plot of LITF verses is achieved. 
Densities are read off directly from this calibration curve 
once LITF is known. The above calibration was used for 
six types of filters viz 

1. Shepp and Logan filter 

2. Parabolic filter and 

3. Four cases of Generalised filter. 

A comparison of calibration curve can be made and 
some interesting points noted. 

Ramachandran filter was also used and found 
unsuitable for the present study. Similar observations 
were made by Kulacki et al (1) in their study of filters. 

Figures (4.1) through (4.6) represent calibration 
for different filters used. Fig. (4.1) is for the Shepp and 
Logan filter which can be compared to that obtained by 
Schlosser et al (2). Fig. (4.2) represents calibration for 
parabolic filter and Fig. (43) to (4.6) for four different 
cases of generalised filter ( sec. 2.4). 
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Table 4.1 Calibration of LITF using Shepp & Logan Filter 


Case i 


1 LITF 

i 

1 ' 
Water • 

1.0 

; 1.075 

j Walnut 1 

0.73 

0.6986 

Pine | 

0.41 

> 0.351 

i 

i Air } 

o 

• 

o 

0.0966 


Table 4.2 Calibration of LITF using Parabolic Filter 


! Case : 

P 

LITF 

; Water j 
| Water I 

1.0 | 

0.751 

s Walnut J 

Sr \ 

0.73 ; 

$ 

0.480 

, l 

■ Pine ; 

o.4 ; 

i 

0.222 

; Air * 

o.o ; 

0.08 


L 


Table 4.3 Calibration of LITF using Generalised Filter 


7™ “ ' 

> Case 

j J 

j Gen 1 

1 

i 

Gen 2 

\ 

L G s n 3 

j Gen 4 

i 

i 

. Water | 

• 

o 

i 

: 0.924 

> 

, 0.915 

« 0.926 

: 0.928 

Walnut | 

> 0.73 

\ 0.736 

‘ 0.728 

- 0.682 

; 0.74 

■ ^ . ! 
Pine 

; 0.41 

• 0.489 

f 0.483 

\ 0.491 

\ 0.492 

V 

3 Air ! 

1 l 

! 0.0 

! o.o 

L . 

; 0.0 

i 

o 

• 

o 

: 0.0 

i ... _ 
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MEASURED VALUES OF DENSITY 


FIG. 4-1 CALIBRATION OF SHEPP AND LOGAN FILTER 
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FIG. 4- 


O AIR 
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2 CALIBRATION OF PARABOLIC ' FILTER 



O AIR 



FIG. 4-3 CALIBRATDN OF GENERALISED FILTER NO 
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FIG. 4-5 CALIBRATION OF GENERALISED FILTER NO.3 
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MEASURED VALUES OF DENSITY 


FiG.4-6 CALIBRATION - OF GENERALISED 
FILTER NO. 4 
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Table 


Table 


t.4 Values of LITF for various scans for Shepp and 
Logan filter 


' v ScajT^fo>J 

0.9 

j 0.8 ■ 

0.7 

j 0.6 

\ 

1 1 \ 
i > 

1.15 

I 0.97 

i 

0.93 

; 0.71 

i 2 ; 

0.97 

i 0.91 i 

0.85 

; 0.74 

3 •: 

1.02 

\ 1-05 J 

0.83 

; 0.90 

! 4 \ 

1.01 

! 0.96 - 

1 i 

0.88 

\ 0.74 

6 ; 

1.23 

i 1.03 : 

0.86 

\ 0.73 


* Rejected 




1.5 Values of LITF for various scans for Parabolic Filter 


1 

2 

3 

4 


ot 

-! 


■ - 

t •' * 

T 

£\ 

0.9 

; o.8 

0.7 

t 0.6 

i 

N 

- m 

i 

3 

* 





1 w 

1 

' 


0.84 

l 0.7 

i 0.69 

j 0.52 
> 

& 

* 

* 

0.69 

i 0.66 

; 0.62 

i 0.55 

* 

i 

/ 

0.73 

$ 

> 0.77 

f 

i 0.61 

; 0.68 

■ { > Rejected 

4 

0.72 

[ 0.7 

: 0.66 

; 0.56 

< 

s 

* 

» 

1 

s 

0.91 

; 0.76 

i 0.63 

j 0.54 

5 

1 



‘ . 

... ; — . *. 



6 
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Table 4.6 Values of LITF for various scans for Generalised 
Filter Number 1 


! 

f * 0.9 

\ 0.8 

i 0.7 

i 0.6 

Scan N 

bxj 

4 

! 

1 

i i 

[ 0.82 

j 0.76 

| 0.684 

; 0.573 

» 2 
! 

: 0.81 

1 

) 0.744 

\ 0.67 

\ 0.57 

; 3 

i 0.813 

' 0.74 

1 0.643 

\ 0.564 

! 4 
» 

■ 0.824 

| 

* 0.75 

0.66 

} 0.543 

1 6 

( 0.84 

i 0.77 

■ 0.694 

- 0.584 


Table 4.7 Values of LITF for various scans for Generalised 
Filter number 2 


f 

Scan 

0.9 J 

0.8 

- y ™ 

1 0.7 

j 0.6 

, % 

' 1 | 

0.81 \ 

0.75 

; 0.68 

■ 0.57 

■ 2 

0.804 \ 

0.74 

‘ 0.66 

| 0.56 

' 3 

0.81 L 

0.74 

; 0.64 

? 0.56 

; 4 

t ; 

0.82 ! 

0.74 

\ 0.65 

| 0.54 

i 

j 6 

0.83 J 

s 

0.76 

■ 0.69 

j 0.58 

_ 

1 


1 

j 
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Table 4.8 Values of LITF for various scans for Generalised 
Filter number 3 



0.9 { 0.8 | 0.7 \ 0.6 

» . * 

0.823 \ 0.76 | 0.69 \ 0.58 

0.82 ? 0.76 ; 0.672 I 0.57 

t ! ! 

0.82 l 0.75 ’ 0.65 5 0.57 

v l i 

0.825 : 0.755 \ 0.663 ? 0.55 

■ I I 

0.84 { 0.775 0.7 \ 0.59 

■ I - -~1 - - - 



Table 4.9 Values of LITF for various scans for Generalised 
Filter number 4 


; Vi 0.9 ; 

Scan__No>^ 

0.8 

‘ 0.7 

t 

z 

‘ 0.6 

* u rv 

1 

: i 

* 

? 0.83 ' 

0.765 

\ 0.69 
> 

: 0.58 

j 2 

r. 0.82 ‘ 

0.753 

: 0.676 

; 0.574 

' 3 

i ' 

- 0.823 ; 

> 

0.75 

; 0.652 

\ 0.572 

' 4 

! 0.833 
> * 

0.76 

! 0.667 

i 0.546 

i 6 

■ 0.843 f 

0.779 

I 0.70 

i 0.59 


t _ j. • 
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Table 4.10 Errors obtained for different scans (Gen 4) 


(Scan No. 

j <f>=0.9 

l<V?=0.8 

t 


j<^=0.6 

s 

i 

V H 

\ 1 

! 0.07 

\ 0.03 

1 0.01 

1 0.02 

i 

2 

| 0.08 

\ 0.05 

l 0.02 

: 0.02 

l 

i 

1 3 

\ 0.08 
* 

j 0.05 

f 0.05 

i 

j 0.03 

i 

i 4 

I 

. 0.07 

i 0.03 

| 0.03 

[ 0.04 

i 

| 

\ 6 

f 

' 0.06 

| 0.02 

o 

• 

o 

i o.oi 

? 

> 

s 

. 

y - „ 

j, - 

- ' w - 


j 


Table 4.11 Relative errors Vs ( Gen 4) 


r <p> 



" i aZ{> 


s 0-9 

l 0.072 

j 0.08 j 

* 0.8 

[ 0.036 

\ 0.045 j 

[0.7 

> 0.023 

? 0.033 } 

; 0.6 

1 0.02 

1 0.030 : 




O SCAN NO- 1 
A SCAN NO- 2 
Q SCAN NO- 3 

1'2 r ▼ SCAN NO. 4 O 

■ SCAN NO- 6 


! i 9 

>- i 

t 1 -Of- 5? A 



</>> 

MEASURED VALUES OF DENS'T V 


FIG. 4-7 COMPARISON BY SHEPP AND LOGAN 
OF ACTUAL DENSITY AND DENSO 
FROM CAT SCAN 
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O SCAN NO 1 
A SCAN NO 2 
Q SCAN NO 3 
▼ SCAN NO 4 
■ SCAN NO 6 


12 r 





<p> 

MEASURED VALUES OF DENSITY 

FIG. 4-9 COMPARISON BY GENERALISED FILTER NO. 1 
OF ACTUAL DENSITY AND DENSITY 
FROM CAT SCAN 


i 



RECONSTRUCTED VALUES OF DENSITY 



0 0 2 0-4 0-6 0-8 1-0 

<p> 

MEASURED VALUES OF DENSITY 


FIG. 4-10 COMPARISON BY GENERALISED FILTER 
NO. 2 OF ACTUAL DENSITY AND DENSITY 
FROM CAT SCAN 




MEASURED VALUES OF DENSITY 

3. 4-11 COMPARISON BY GENERALISED FILTER NO- 3 
OF ACTUAL DENSITY AND DENSITY FROM 
CAT SCAN 




MEASURED VALUES OF DENSITY 

FIG. 4*12 COMPARISON BY GENERALISED FILTER NO- 4 
OF ACTUAL DENSITY AND DENSITY FROM 
CAT SCAN 



Tables ( 4.4 to 4i9) show the reconstructed 


density values for various scans and for various filters. 

Fig, (4.7) to (4.12) includes a comparision of the cross 
sectionally averaged values of density by the scan method 
to other method, denoted by/^2 g and ^'^/respectively. 

The darkline indicates a perfect match between the values 
of density measured by the two different methods. A 
noticeable feature in Fig. (4.7) and (4.8) is that the 
values of s are greater than the values of<^/v>for 

all the scans. However, for all the cases of the generalised 
filter (Fig. 4.9 to 4.12) we see a considerable improvement in 
in the LITF values compared to those of the Shepp and 
Logan, and the parabolic filters. 


Table (4.10) summarize the deviations of <^/^. 
values from values. Fig. (4.13) represent Table 

(4.10) in a graphical form. There is an under estimation 
of true values by the reconstructed density values as can 
be clearly made out in Figs. (4.9 to 4.12) . The difference 
S</> , between these two values follows a pattern with 
width of /\,< ; ^>value. 



The average of the absolute deviations of the 

j v i 

the values of < ^ i.e> the average of j fK^/^jfor various 
scans is shown in Table (4.11) and graphically represented 


in Pig. (4.14). Also appearing are the relative errors* 
4^^ f/ /Which are plotted against P in Fig. (4.14). The 

~7'}\ ; 

^ -I / ~ 

quantity IJ is computed by the following equation. 


<: \ jf- sp\{ f 

> °c ^ / I/-’, 


— v • i / n 

where n is the total number of scans considered. Since 


the error band of values remain fairly constant for 

v •- s 

all density runs, a decrease in the relative error with 
decreasing density is observed. 
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MEASURED VALUES OF DENSITY 


FIG. 4-13 ERROR IN DENSITY MEASUREMENTS 
FOR GENERALISED FILTER NO. 4 
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MEASURED VALUES OF DENSITY 


G.4-14 RELATIVE ERROR IN DENSITY MEASUREMENTS 
FOR GENERALISED FILTER NO. 4 
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CHAPTER— 5 

CONCLUSIONS AND RECOMMENDATIONS 
5.1 Conclusions: 

The following conclusions have been drawn based 
upon the results of the present study. 

Meaningful results have been obtained using the 
techniques of computerised axial tomography (CAT) for density 
measurement in bubbly air -water flow. For cases investigated 
in the range 0*fe^j^l.0 , the uncertainity band in the 
measured values of average density is found to be -0.03 , 
with respect to density determined experimently by direct 
measurement. The density is underestimated for all scans. 

The dotted line in Fig (4.9 to 4.12 )<(/) =1 .07 cl early shows 

that the error can be overcome by a general approximation. 
This is true for all Generalised filters. 

The results of density (scan) values using Shepp 
and Logan filter and Parabolic filter are enclosed to give 
a proper insight into the problem. The improvement the 
generalised filter offers can be high lighted by viewing 
Figures (4.iy to (4.12). 



The above method is approximate and its 
limitations can be exposed by using less number of scans (N) 
and less number of rays per view (LITM), Proper cace should 
be taken while choosing N and (LITM. The values of N=20 
and LITM = 40 was used for the present study. 

The method ( 'New Algorithm'. ) is accurate and exact 
but takes longer time to evaluate. The method also removes 
the necessity of rotating the source provided the distribution 
is centrally or. radially symmetric. However repeated scans 
at different scan angles will be required if the flow under 
investigation happen to be radially assymetric. 

5.2 Recommendation: 

1. The Generalised filter used needs further investigation 
and selection of various values for Pand Zeta. The 
mathematics part of Generalised filter needs proper 
understanding. 

2. The Noise error analysis as suggested by Kwoh et al ( ^ ) 

is to be considered and this is an important factor. 

3. Plexi glass was used in the set up for data calculation 
but its effect was neglected i.e., not taken into 
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account in the algorithm. Suitable modification 
in the algorithm is necessiated and 'New Algorithm 1 
comes into prominence for this problem. 

4. The 'New Algorithm' explained in section 6.0 

has to be adapted for experimental data by suitable 
modifications . 
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CHAPTER - 6 

A NEW ALGORI THM 

The Algorithm described and used in the present 
study is based on convolution i.e., filtered back projection. 

It is used primarily because of the high processing speed and 
reasonable engineering accuracy. 

A " New Algorithm " is being proposed for the reconstruction 
studies in the two— phase flow area. To test the concept we 
generated data from a known density distribution /then used 
the general data to reconstruct the assumed density distribution., 
A brief summary of the procedure is given below: 

Step 1 Data is represented as an . array and given as a 

column of n elements where, each element is the 
chordal averaged value at different scan angles. 





We note that N or LITN explained before and 
defind as Number of scans is unnecessary in this 
case. 

Step 2 s~ We assume a function g, a polynomial of degree 
'\\7' n in r 2 . 

2 2m-2 

where g=x 1 +x 2 r +...+x m r 



B COS& 


<x 2 +r 2 sis 2 &> k 



Pig. 6.1 Assumed distribution for New Algorithm. 
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where x is the coefficients of powers of r 2 . 

x can be represented as /x, \ 

/ 1 ' 

I x 0 

i 2 . 

I 

X= S . 

I . 

\ X -n 

\ m 


(6- 1 - 


Step 3 Nor JX = D 

where J is the weight of each part which when multiplied 
by x will give rise to d. 


where 


(i) 


j _ j 

= 2R 2 Cos 9, ! R 2k “ 3 + (k-1) Sin 9; 

L J 


2 k-1 


k= 2(1) m 
i= 1(1) n 


(i) 


= 2R CosQ. -th ^ _ - . 

i , x element or 1st column 


i = 1(1) n, This gives n x m dimension for J» 


for derivation of J see Appendix E. 

On solving the system of linear equations J x = D, we obtain 
the X-vector 

Step 4 5 - The reconstructed value at the point r is given by 

hi 

F(r) = ^ T X = r2lx i 


where 


P ■/ ? 


^ _ 
it 0 


*2m-2 
r 


Step 5 s- The data used is generated especially to verify the 
'New Algorithm 1 . The data element d is the chordal average 

Y&jfi C7 

and is calculated using the formula d - 2 j f(r) dx. 
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where £(r) — f ( ; x^ + sin%" ) (from figure 6.1) 

H 

The accompanying data provide ample proof for the 
success of the elgoirithm* 

A FORTRAN program has also been written and employed to 
give the best of results • NAG library routine has been used 
to perform Numerical Integration in calculating data and for 
solving system of linear equation to find X by Eqn (6,3). 

Several functions of f(r) were used to test the algorithm. 
Some typical cases are : 

f(r) =;"1 
l r 

J exp(r) 
j exp ( -r ) 

Tables 6.1 to 6.4 show the reconstructed values, F(r), 
and the exact values of the various function in the interval 
0.0^ x 1 .0 in steps of 0.1. 

We observe that the n ew algorithm predicts the unknown 
fuction very accurately if it is radially symmetric. This algo- 
rithm has to be modified (a different J-Matrix has to be used) 
for use in radiation attenuation methods. Also, the algorithm 
can be extended for sitation when the fuction is not radially 
symmetric. Such a case will involve multiple projections as 
currently being practised in CT algorithms. 
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TABLE 6.1 



.... f(r) = 

= 1.0 _ i 


jsi No. 

* 

.! f (r ) | 

^ i 

; F(r) ! 

d j 

i 

■1 

; l ,ooo 

1.000 5 2.000} 

2 

\ 

! 1 .000 ! 

; 1 .000 1 

1.978 j 

: 3 

i l .ooo 

* i 

j 0.999 | 

1 .912 j 

.4 

1.000 

0.999 : 

J 

1.805! 

! ! 

1.658j 

• j 

•5 

; l.ooo 

j • i 

1 .000 

.6 

l 1.000 i 

j 

i 1.000 

1 .475 i 

(7 

t 

1.000 

t i 

; 1.000 1 

1.259f 

■8 

1 1.000 

| 1.000 | 

1 .015 j 

• 9 

’ 1.000 

: 0.999 j 

0.749? 

} 

: 10 

j 1 .000 ; 

! i .ooo 

; i 

0.467? 

: : ll 

j 

j 1 .000 i 

I 0.999 

0.174 \ 

TABLE 6.2 




fir ) 

= r 




Is 1 

No. ' 

fir) 

l 

F(r ) ' 


d 


1 

i 

j 

] 

0.000 

1 

0.095 

i 

if 

J 

1.000 

i 

| 2 

i 

0.100 

< 

0.119 

l 

J 

1.045 

f 

1 3 

! 

0.200 

> 

0.186 

t 

t 

1.118 

t . 

f 

-4 


0.300 


0.284 

t 

\ 

l 

1.178 

t 

-5 

] 

0.400 


0.398 

? 

4 

1.200 

5 

\6 

t 

! 

0.500 

J 

0.508 

* 

t 

1.170 

f 

|7 

| 

0.600 

l 

if, 

0.605 

i 

t 

$ 

5 

1.076 

l 

i 

|8 


0.700 

j 

0 .695 

j 

0.923 


-9 


0.800 

i 

0.798 

i 

0.713 

* 

•10 

f 

t 

0.900 

■.J. 

0.905 

1 

j 

0.458 

5 

j 11 

l 

1.000 

\ 

1.017 

1 : 

0.174 

l 



TABLE 6.3 


r _____ 


fir) = 

exp(r) 



! SI No* 


fir) 

i 

P(r) 

T " 

d ! 

f 

: i 

: 

1 .ooo 

\ 

1 .095 

i 

J 

3.436 J 

;■ 2 

1 

1.105 

1 

l 

1.124 

j 

3.473 | 

’ 3 

J 

1.221 

f 

1.207 

! 

3.518 J 

; 4 

* 

J 

1.345 


1.334 


3.522 j 

: 5 

f 


1.491 

< 

1.489 

^ - 

3.450 | 

! 6 
t 


1.649 

* 

1.657 


3.273 j 

| 7 


1.822 

* 

f 

1.827 


2.967 j 

i 8 


2.014 

* 

2.008 

1 

2.521 : 

: 9 

t 

j 

2.225 

‘ 

2.224 

i 

| 

1.941 ■ 

: io 

f 

2.460 

1 

2.465 


1.246 | 

1 li 

( 

Jr ^ 

2.718 

i 

2.735 


0.473 | 
- , J 



TABLE 

6.4 





f (r) 

= exp(-r) 



Psi No. 

► 

1 

fjrT ^ 

■ — 




-4- 

" ~ \ 

■ i 

\ 

i 

1.000 


0.967 

! 

1.264 5 

§ 

i 2 

t 

i 

0.905 


0.915 

i • 

1.205 ; 

| 

. 3 


0.819 


0.816 

& 

* 

l 

1.091 j 

J 4 

X 

0.741 


0.739 

1 

! 

0.944 | 

: 5 

j 

0.670 


0.673 

i 

| 

0.811 \ 

? 6 

i 

i 

0.607 

S 

0.605 

i 

0.671 ; 

• 7 

■ ■: ■ 

0.549 

i 

0.549 

> 

i 

0.536 j 

; s 

? 

0.497 


0.497 

* 

* 

0.409 \ 

| 

( Q 

5 

t 

0.449 

- 

0.449 

v‘ 

5 

s 

0.289 j 

i 10 

i 

t 

0.406 


0 .406 

J 

j 

0.175 j 

’ 11 

( 

t 

0.368 


0.368 

* 

0.064 j 



T’nis 
f. an 

to o 
i:'ITp; 


= d ? ta t0 the tomography system it accepts 

2® 9 ?, „?t? foT na rallel to convert w.hi rh converts the data 
arailel beam. 

GER M , LITN 

WMMAX,DFLCiAM,T,BRTA,H f LAMBDA, THETA 
o l 497 i 4 tance from soarce t0 center of object 
the diameter of the object 

is the number of projections in scantfan beam) 

~20 

the number of rays per viewCfan oean) 


tr y J. a tes maximum lambda 

LAMMAp r/(2.0*SQRT(1.0-CT/(2.0*0))**2.>) 

spacing in terms of iamoda 
u a, LL AMs2 . U*LAMMAX/FLQAT(M ) 
valculates angular spacing between views 
ftLP,(4=3. 1415927 /FLOAT t ABSiL Il’N) ) 

the output in f 1) e named input. iat 

WKITE52l")D!^ FILfc:a ' ' TNt, ‘J T '» A ^K3 Ss 'R^‘ JT 'i 

I'.’dI TE ( 2lr*) LAMM AX , DELI. AM, ALPHA 
WRTTE(21 ,*) M# LITh 
CLOSE (UM I i’=21) • 

Tf PE 300 

FDRMATC' NEXT PROGRAM is C3LV£RT f ) 

SFOR , * 

EMD 


maximum 

1 r»*StiP’ri 



Uflti lA Cl 1 1 tl lit* 


7 


.riis program converts fan- beam data into parallel beam 
Jata. tne data are then passed to tne filter function 
JtoTta.T' (oicn is geometry dependent) 

I id T f* flL.lJ L 4 !•. f 'f * iVa . T f mil. m* • y t « » 


I : i T h G h R si . b X I' n f L T T M , f ! . K , J . R , I 
"k r '\' r I'SlJ '■ l f IiAMMAX , LIT A , f)ELL 


x LA*MjITH( 41 ,41 ) ,T 

K m A r ■ I LI IPU 1 . 41 ) ,BFTA,LAMBDA,LI 7 L,THETA, 3 A A K A 
i c - n" , II? ? I I s ? | f F Il.K* ' I W p UT . DA r ' , AC C ESS= ' SE 2 1 i ' ) 
F jRMAT( M S . 8 ) 


fnc data is read from file inout.dat created oy Input. for 

RE AD (21,*) D, T 

RhW\!»(2i ,*) bAMMAX,Dfc:T.LAM, ALPHA 
RdAuf 21 ,*) M.LITN . 

CLOSE ( IIM I T=21 ) 

M Is the number of prolections in scan (parallel) 

11 = 2 0/. . 

L l T M Is the n timber of ravs per prolection (parallel) 

L I 1’ M — 4. 0 

fx gives the option to specify the required data file on 
The tty up to 5 letters maximum, e .a* water , wanut , etc . , 

ft it A D ( 5 , 9 ) FM 
FORMATE 2A5 ) 

00 10 J = 1 , LITI4 + 1 

OP CO CUNT T=5 3 , ACCESS* *«EQI M * , F I LEsFN ) 

READ (53,5) CLITH C I , J ) , 1 = 1 , M+ 1 ) 

CLOSE (UNIT=53) 

CONTINUE 

LITU is the fan beam data as an array 
DO 3 00 J*1 ,LITN + 1 

Tne next three statements change LI T H to density values 
00 100 1=1 f LITN+l 
LETHCX, J)sf .0-LITHCI#J) 


U U X V V X l I i j U n T * 

LETHCX. J)=f . 0 -LITHCI#J) 

CONTINUE 


LEVA is the distance between parallel rays in each view 
LITA=T/FLL)ATU,TTM) 

CONVERT TO PARALLEL 

00 80 K=1 , LlTHtl 

L i T f , = L I T A * F 1, 0 A T c K - ! ) - c T / 2 . ) 

GAM 44= AS EU ( LITL/D) 

UAMBOAsOtC SIN (GAMMA) /COS (GAMMA) ) 

S=1 ^M/2 1 1 FIX (FLOOR (LAMBDA /DELL AM) ) 

IF (S.RQ.O) S=1 
D f 3 7 0 J-i Nfl 

THETA* ( 3. i 41 5927/FLOATC N5 )*FL0A’F(J-1) 

BETA=T»ETA+GAMMA , aj4 _ n __ 

I=IFIX (FLD5r(BETA/ALPHA))+1 

LI TP is the parallel beam data as an array 

Interpolation is * being done in the next few statements 
IF^ (S^LE. M)LITP(K?J)=LITP(Kr J) t (LAflBOA/OELLAM* 

1 FLOOR(LAMBDA/DELLAM))*(LI|HCS+l,I)-LirH(S,I)) 

IF (I.LE.L1TN) LITP(K, J) = LITP(K|J) + (8ETA/ALPHA-{'L0DR 
1 (BETA/ ALPHA) }*(LITHKS,I + i )-LITH(S, I) ) 


CONTINUE 

CONTINUE 




Hatton in file named qeoro.dat 

<RI TE ( 23 U) iJlA 
?QRHAT<Fl6.9) 

:i.OSS (UN1T=23) 

i>£j;te parallel data in file narad convrt.dat 
3 P k- N „ ( U N I Tb 2 4 , F I l,e* ' C ON V R f ' , A C Z E S S = ' S E OOUT ' ) 

.RITE (24,6) fj.LITM 
FORMAT(I2,3X,l5) 

V< ITE (24,7) ( (LITP( J , I ) , J=1 , LIT 4 + 1), 1 = 1, 7+1) 
P0RMATC4(E15,8.3X)) 

CLOSE (UUIT=24) 

TfPE 410 

FORMAT ( '-NEXT PROCRAM IS FILTER' ) 

ST Dp 
SOD 

Tots function coitioutes the greatest Integer less than or 
euual to the argument, it returns a- real result, 

FUNCTION Fl.DORCX) 

REAL X 

FLOORalFIX(X) 

IF ( FLOOR. GT.X) FLOQR=FLO0R-l . 0 

RETURN . 

END 



/ 'rt, 1 r4l)*LITK<41),8lH f UTC(41,4n 

TdThaER LI Pto,N, J,K,L, 
p jR*tAT( 12 . 2A. 1 4 ) 

^gaa Parallel data from convrt.dat created dv 
OP Erl ( II *T£*?A , FILEa "COHORT ' , ACC£SS= f S E D X ^ ' ) 

R t-J i\P (2 1,*) 'J.LITM 
FOR »! AT ( 1 2 , 3X , 12 ) 

READ (24,*) ( Uj 1TP(K, J) ,K=1 .LITM + l) ,0=1 , Ul ) 


cdnvrt. f or 


F i.) R »! A T ( 1 2 , 3 X , 1 2 ) 
r.o cir r ii.itt-OaS 


f U. HlU 3 , 3 i JAJ J 

CLOSE (JrIlT='j4) 

Read t curler filter from filter.dat created 
OPSd ( U 3 11 = 5.1 .FILLs 'FILTER * , ACCESS® * SEPl d ' ) 
L X JP 1 is tae filter function data 
0 1 30 K = 1 , MTM+i 
R £ A 0 (St.*) HTwi(K) 

F 0 R A A T ( G i 5 , 8 ) - 

F iniAr(4KiS,«) 
c-Q.vrii'HJE 
CLOSE (UMIT=51) 

Compute convolution ■ 

03 70 L = 1 , LITM+1 
03 60 J=!,M+1 
SHMsO.O 

0-Q Sa nK*1 , LITM+1 

SUMsSUMf biTPCK, J)*LITW( ABSCL-KJ+l ) 

CONTINUE 
LITCCL, J)sRUM 
C 3 N 1 1 WHS 

CONTINUE , „ 

Write convoluted data in convlv.dat 
OPS) ( U M IT a 5 2 , F l h E = 7 C 0 N V T , V ' , A C C S S S = 


by filter. for 


K) + l) 


OPE;) ( LKHT=52 fFTliE* CONVliV , 
writing tne convoluted data 
WRITE (52 , * ) M . 01 TM 

f%. <«* 4 'A It mm. 4 1,1 1 4 


ACCeSS='SEOQtJT f ) 


1)3 100 KaJ.fJ + 1 
tRITE (5 2 , * ) ( LITC ( ti , K ) , L*1 
C OUT Id US 
CLOSE ( !J N I T* 52) 

WE 200 

FORMAT! '-NEXT PROGRAM IS 
STOP 
Eli) 


1 ,HTM + 1 ) 

, "SUPERI’4P0SITnw H 





U S T f L Y T v C ?U- l i rLITC(41 ' 41>/SHM # r, THETA, 
l 4 l M:itn J , U t X | l . bl T K , N 

' 7 data from coavlv.dat created by convlv.f or 

l®n Jsi JrSlLlI?*' '™*W'KCESS*>SeQI»') 

FORMAT(I2,2X,I2) 

LU'C is the convoluted data 

■) J b (i=l,N + l 

CONTINUE 
CLOSE (U 011=54) 

calculate super imposition of data 

I ) . ..1 J G X. ’-"I# ,u I f M t % ' - 

DO 20 Y = 1 , tITf'Ht 

SUMsd.O 

DO 10 J = l ,N 

THETAs (3, H 15927 )tFLOAT(J-n/FLaATC.’l) • 

TRUET=FLOATC X-LITm/2-1 ) ♦COS ( THET A ) +FLQ AT ( Y-L1TM/2** 1 ) 

1 ♦SIN (THETA) 

T = TRUET + FLOaT ( LITM/2+1 ) # 

L=IFIX(FLOQR(T)) * 

SUMsSUMf ( FLOQRCT + l „0)"*T) *LtTC (L , J ) + (T-FLOOR ( T) ) *L ITC ( l+l » J) 
CONTINUE 

L ITF Is the bacK-pro jected data as an array 
TjITFCX. y)*SUM/FLOAT(2*N) 

CONTINUE • 

■ CONTINUE 

('(rite output data In supimp.dat 
LIT M*LI TM + J 

OPEN (UNIT=55fFILE='SUPIMp',ACC6SS-'SEQ0Ur') 

WRITE (55.2) LITM.tlTM 

N K I T £ ( 5 5 »1 ) (CLITF(X,Y) ,X=1 ,LITM) ,Y=1 ,L1TM) . 

FORMAT ( 4E1 5 , 8) 

CLOSE (UNIT=55) 

TYPE 1000 

1 00 FORMATC '-NEXT PROGRAM IS "PRINT"!') 

STOP 

SMO 

FUNCTION FLOORCA) , , J . . . 

Function is similar to the one explained oefore in canvrt.foi 

R2AJ, A 

JTF 0 (FLOOH^t!a) FLOOR*FLOOR-I.O 

RETURN 
END 



u J Kf I L n h LvSN 

REAL LI TFC 41, 41) 

L , >t XMAX, YMAX,SUPERI, PRINT, PUT, TITLE (20! 

/ i1-t-?-lc 0 ST?13 p l mp » dat created by suoinp.for 

mao' ess {j 5 5 t!iii£?® PI ** ••‘:cES5=- S Egi.-) 

FRRHftTCl2.2X.I2) 

CLOSE CUfiITr55) 

Print data 
Y=9999,C 
X=-9999.0 
DO 30 L = 1 , XMAX 
DO 20 J s t , Y H A X 

IF C X . G’C . LX TF CL , J ) j Y=LI TP CL , JO 

Z=LirF(L,J5 ' 

continue 

-j; jf Pr 35 

FORMAT! '-NEGATIVE OF IMAGE(0=NO, 1=YES ) * 0 
ACCEPT * , ft NTS 
TYPE 36 ■ 

FORMAT ( 'PRINT IMAGE ON TF.RMIN#\L?CO=NO, 1 = YES) ' ) 

ACCEPT * , PRT ‘ ' 

Open print file 
EVENS. FALSE. 

OPHM (U H'fsSb , FILEs 'PESULT ACCESS® ' SEQOUT ) 

00 55 LSI, XMAX 
no 40 0=1. XMAX 

liltiKC J)*IFIX CO. 5 + 99.0* C LITE CL, J)~Y)/(Z-Y) ) 

I F {ft MS. ME. 1 ) LINECJ)=99-T.INE(J) . 

continue 

WRI TE C 56 , 45 ) CLINE(M) , 4=1 , XMAX) 

FORMAT ( 1 H ,5012} 

IF C EVEN) GOTO 300 


WRITE (5 6, 200). (LIKE CM) ,M=1 ,YMAX, 2) 
WRITE (5 6. 200) CLINE(M) ,YMAX,2) 
F0RMAT(IH+,25CI2,2X) ) 

cotq 500 

V WR ITEC 56 , 400 ) (LINE(M) ,M=2,YMAX,2) 
WRITEC56, 400J CLlNF(M),Ms2,YMAX,2) 

n n -iiisdtf i u .a- tl f v T O "4 1 


n n Ji .4 ifcHi \ i# \j f i v wj % r. f \ w / § ri — 

FORMAT! I H+/25 C2X, 12) ) 

EVEN®. NOT. EVEN 
IF CPpf.EQ.l) TYPE 45, CLINE 
CONTINUE 

forma£c/ |6 m1ni&um LITF= ' , El 5 
CLOSE CUNITs56) 

STOP 

END. 


(LtNE(M) ,M = 1 , YMAX) 


MAXIMUM LtTF=* ,E15«8) 




? r3 iit a ^r>H n £?? rates ttle values or Inf ed 

,*itl ? lves a number for t.ne 3ix M valu 
FTFR r NM“41 be cnanaed by using a nigner valu 

r itfir'r2S? n L ,|FH ffi . result.dat created oy priat2 
i V J £ { = 2 1 t\ lLfcs * RESULT ' , 4CCESS= ' SEP I M ' ) 

21 ,*) C CB(l,d) ,Jal,NN),I = l»N»i) 




/It m r'n?»r a *nH2l* u ^ a l es S 116 Courier filter function. 
^ r '- 1 p , J.} ! ^ r „t'unction is the parauollc filter finction 

; l}, y „ }i l 1 A i ~ Ul N 5 1 . W ( 4 1 ) 

: T ^«*-R MTH,M,GEOM, FILTER 
Mad geometry factors 

!lnn }5iJi)? i !t® fc *' eM " ,, ' w::E8s *' se0l "'> 

MRM/1 1* ( 12 , 3X# 12) 

MAO (21. *) Lit A 
’JIM AT (El 5. 8) 

:lOSEU!NIT=23) 

Calculates constant 
H 1, ) = ( 3 . 1 4 1 5927 ) / 1 3 . *LI T A*LI TA ) 

•Pi £ T ~Z'\ . 0 / ( 3 . 1 4 1 59$7*LITA*LIT A ) 
irite filter values 

IPRN ( UNITbSJ , FILE® 'FILTER * , ACCESS* ' SEOOUT ' ) 
fRITG (51,10) (W(L) ,L«1 ,41 ) 




mmmm- - ^ 


Tfl# Lf ,L for ^eneralisei Filter for no 

I'}® ?Mh .nnfiL^ ta ^' 1 p are specified in tne Program 

is the HkG Horary routine DOiAKF 


Tne routine is ca 

INTEGER 14(600) 
r>T '-i f-; m s i □ i-j 4 (^ 000 ) 
E X f i:' K is A L F .1 , F 4 
C ).«mon :s,p,zeta,p 

A = 0.G 

P(F: = 22./7, 

R=2.*PIE 

£PSA8Sa0,lE-0i 

EPSRELsO.lE-Oi 

It A J u=0 

ZETAs4.75E-03 

P-2.0 

DO 10 0=3,41,2 


caned twice for tne two separate term 
0 ) , RESULT (42 ) , VALUE (42),FI*AL(42) 

,PIE 


■■ ■ J A V# f I 4, 

Tne routine below does the first part of integration 
this uses a function called f.3 which is external 
CALL DOt AKF(F3,A,h,EPSABS,F;PSREL,RESULTCN) , ABSERR, 

~L^?,S' )0 ' IW ' 6oS 't FA i L > 

CONTINUE 
on 11 11 = 3 # 4i 2 

Tne routine feelow does the second part of integration 
tnis uses a function called f 4 which is external 
CALL DO 1 AKF( F4 , A, PIE.EPSASS, EPSREL » VALOEf 4) , ASSERR, 

1 14,600, iFAlt.) 


,2000 


• n , £ « V V , A w , o v u , A r H i ll ) 

F 3 R M A T ( E 1 5 • 6 ) 

F 1 N.AL ( M ) a RESULT (N) +VALUECN ) 
FT0AL(N)=FfNAL(N)/(PIE*0»1625*0.l625*N*N) 

WRITE? 2 4, 2) (FINAL (I), 1 = 3, 42, 2) 

4RITEC5,*) I FAIL , ABSERR , A,B 

STOP 

EMD 

FUNCTION F 3 ( X ) 

This is used for the first part of integration 
C'lMFOL N, P,ZETA f PIE 
SUM t = 0 , 0 
AS=0 . 1 §25 
K=(N-3)/2 
03 100 J=i,K 
Y=X+2.*PIE*FLQAT(J) 

Y SuW^S0Ml+YSEXP(-zlTA*¥i) 

F3*SUM1*C0S(X) 

RETURN 

END 

FUNCTION F 4 ( X ) 

This is used for the second part of integration 
COMMON M, P , ZETA , PIE • ' 

SU*-l?a 0,0 ’ 

YY=X+PIR*(FL0AT(N)-1.) 

SUM2=COs[x)*YY^EX^(~ZETA*YYl) 

F 4 =SUM 2 

RETURN 

END 



uaifilc are tr >e Generalised Filter tor t*:ven w's 

If,} r^?,*? f ^ Zeta J a ? d p are specified in tne Prouram 

?Hr I^M% used the nag liorary routine oofAhF 
^r,h tWtSOO) 

; |Sipn 4(?DO0) , PKSULTC42) 

H N A L lr 2 

Ou id , P , ZETA, PIE 

11 * . .. .. v ' 

the lower limit and B is the upper limit of the integral 


and fpsreL are the absolute and relative errors 


f FAIL=h 

ZST A = 4 « 75R*0 3 

Zeta is the damping coefficient 
0~J U N = 2,40.2 

CAIjb f>Ul AKF(F2,A,B,fe;PSAB5,EPSREl,f RESIH.T( N) , ASSER 
1, d,2i)U0,IW.600, IFAIL) 

R E s U h T ( N ) s R E S U h T ( N ) /• c p I E ■ * 0 , l 6 2 5 * 0 . 1 b 7 5 * N * N ) 

CONTINUE 

Wf<ITE(24, 2) C RESULT C I ) , 1=2 ,40,2) 

Result nas the value of the integral 

F jRd/vX ( El 5 .6 ) 

SIMP 

FUNCTION F 2 ( X } 

These statements are the value of the integral 
employed In this Program,!#. r .t. x 
C Mviofi M,P, ZETA, PIE 


162S 


uu i u j - i t r\ 

YaX+2.*PIF»F00AT(J) 
YI=(Y/(AS*ZN))**P 
SUM = SUM+U*EXP(- 2 ETA*Y 15 ) 

F2=SUM*C3S(X) 

RETURN 

END 


>•& 



rnese are 
Tme value 


X + FXPC 





t- : ip S vaf SL nf generalised filter for n=l 

E 1 ,}® ^w^nf? £ 2eta P are soeclfled Ln tne prog 

T^f PGER r I^cloo ) US€d the aaiJ ^Prary routine dot 

RSAl”w(2000) 

EXTERNAL FI 

k~‘).D ' , 

3 - 2 3 * / 7 . 

E P S A R S- 0 « 1 E- D v 4 
E?SREL=0. IE-02 
t F A I L = u 

CAU. HOI AKF(F1 , A, B,EPSABS,EPSREL, RESULT, A35SRR, 
i. *1.2009, 1 W. 600. T FAIL) 


h 6 s u b f 


'Mbmmaio. 


1625 ) 


WRITE (5,*) RESULT , ABSBRP , A , B , IF A JO 

STOP 
F. I 0 

FUNCTION FI (X) 

Toes e statements calculate tne value of tie Integr 
function used and says than is w.r.t. x 
ZSTA=4.75E-03 
p-2.0' 

Ft =X*COS(X)*EXPC-ZETA*(X/G.l&25)**P) 

RETURN 





a e s e f 
Odd an. 
altera* 
in K o r ‘ 

0PE«(U1 

OPEK(Ot 


vSn J 'iH a m,tF‘L val ^ s 01 Generalised filter's 
V i?npOfS», S f^ e ? ir * a sin, 3le fiie inserting 
air inputs. Tne output is available 

; -7i trrpl2EVL^r,9 u g n Jdd.dat and Sven.dat 
,Zi\ ArrPcf^ -flQi II r £11^- 'EVEU ') 

■’^2^ACChSS='SEQ I N # , FILE= 'ODD', 

k 



Scan No. 1 - Normalized 















Table B. 8 Scan No; 2 - Normalized 
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Tnest are tfie values for Scan 6 of 701 ( 0 * 0 , 7 ) density for Generalised filter 
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This Program Generates the Data tnat Is to bo employed 
in the 'New Algorithm' for testing' it. This uses a WAG 
library routine DOlftKF which does Integration over any 
defined limits. The function is to be specified oy the i 
EXTERNAL* F • 

DIMENSION H t 8003, XHC 200 ), KES0LTC1 1 ) 

COMMON TH(ll) 

sgEBisa *s?? 4 *assf?f =;3«"' ' ri, * s “' ruBT ** * 

A. = 0.0 

DO JO 1=1.11 
B = COSD(fHhh 
EPSRGLasO. 00001 
£PS&RSs50 .00001 

(?Af.L D01 $KF(Fl A, B,EPSARS,6PSREL, RESULT!!)# 
ABSERR#W,800,t»N,200,IFAIL) 

CONTINUE 

DO 29 1=1,11 

RESULT ( I )=2 .OtRES'UjT Cl ) 

OPEN C UNIT*2 1, ACCESS* 'SEQOUT'. FILE* 'DATA') 

MR I TB£2 1,9) C RESULT ( I >,1*1,11) 

FORMAT C llEli ,5) 

STOP 

END 

FUNCTION FIX) 

COMMON' THC J 1 ) 

XffWIlfStrJl' ' 


F2f§0Rf(XfX+t*y) 

F*EXP(*F2) 


RETURN 

END • ' 




on a a 


9 

10 
50 


c 


This Program Generates the Data for the 'wew Algorithm' 
The MAG library routine used is the noi AHF routine. 

Tne Function is to oe specified oy the user. 

PROGRAM FUR INTEGRAL EQN 
REAL A , B, EP5R, Ft£LERR,TH , DT 

integer ifaii.,nlimit,n 
REAL D01AHF 
EXTERNAL FUN 
COMMON TH 
A a 0.0 
NLXMIx a 0 . 

EPSR s 1.0E-5 
IFAIL f 0 
DTs 8.5 . 

oo io i 3i,n 

TH a FLOATCJt-i ) *di 
R ' s C* 0 0 C Y §■§ 

ANS a D01AHF ( A , R , EPSR, N , RELERR, FUN , NLIMIT, IFAIL ) 

ANSr2.*ANS 

WRITE C2i *9) ANS 

FORMATCHEil.5) 

CONTINUE 
STOP 

FORMAT' (5X, 'THETA*' F4. 
i 'ESTIMATED ERROR* 'E1S, 

END 

REAL FUNCTION FUNCX) 

REAL X,K 

COMMON TH • , 

K a SIND(TH) 

FUN a EXP(-(SURT( (X*X) 
fUiial.O 
RETURN 
END 


1, 5X,' INTEGRAL*', FB, 6, 10X, 
5,3X,'IFAIL»<,F5.3) 


+ tK*K) ) ) ) 


t»V'- 

f. 




DUO 


mis Program Generates Data for specific t(r )=r in 
order to test the 'Mew Algorithm' ,'Tats doesnot utilise 
the MAG library routine . 
dimension th(Il) f r(U) 

open(units21#access='segin',flle*'tneta') 
read(2i,*) (thC i) , i=l , 1 1 ) 
do 10 isi ,11 

Y*Cl+cosd(th(l))}/sind(th(i)) 

10 r Ci)=co5d(th(l) )+( (sind(th(i)}*sin1(th(i)))*alog(y)) 

open(uhits23,access«'seqout ' , £ ile='newdat* ) 
write(23,9) {rCi),i=l,il) 

9 f ormat^ lien ,5) 

stop 
end 



c 

c 

c 

c 


IUnV! r ?h?!® solves for x using tne relation jx= 0 
Tnf is calculated using tne formula for J 
ais0 * nc ludes printing of f (r) l,e.. 
Function Values and errors in Reconstruction E. 

D 
D 

LOGICAL SVD 
N»U 




X { 1 1 ) 



WHDC12-J) 


■ - 

MB 


SESWHP 1 !*^ access* 's 
^ltiiHMl S! li4SCEssa;ssgiM;' f FiLe* , poi«Ts 
OPEfi i !JNIT»23# ACCESS* SEQIN ,FILE* f DATA'} 


N'f FILEe'THETA'} 

^ * PTr.ffifDrttAlfC* r 


READ (2 
READC22 * 
READ (23 i* 


DO 10 J>i,U 

WOT” 

00 JO' Jtt2 N 

(FLOAtCjp.l, ))#(!.+ 

1 Ly-i) * ( i . -o*o )/q) *c ( j-i , i ) ) ) 


T0Ls*S.*10**C~4) 

00 19 K *1 ,11 

dS K 25 X i 2 l 

00 29 J* I , M 
A(t, J)*C(J,I) 

OPEN jUJN IT-2 3 ^ ACCESS* SEQOUT' , FILE** RESULT ) 

Ips0 2 ' 9) 

CALL F04JGF(ll,il,A, 11,D,T0L,SVD, SIGMA, IRAWK,W,44,m 
NRITEC5,*) SVD, D, IRANK 


iTed Reconstruction output 


OPEN 

WRITE(55 f / / 
FORMAT! 1 lEJ. * w 

FORMATC 1 E 1 l.SL 
CLOSE (UNIT* 55 ) 


ESSs'SEQOUT' } 


STOP 

END 




p E?2 r ?™_ chec * :s 'the acc 

nff'^crAv H«£? r ST routine n 
DIME^SIOh- hi 1 1) , B { 1 1 1 1 1 ) , c { 

gP|^U ? »XT=23 f ACCF, s Ss*SRQi.\ ! * 
RfiA0(23,9) CC8fI,J), I=l,li) 
READ(23,9) (C(I) ! 1=1 ,11) 

DO 10 1=1,11 
DO 10 J=1 ,11 

MIU72l?9]Y ccrs *‘'* fm,T 

FORMAT ( 11 L‘l 1 ,5) 

STOP * 

END 


FILE 





nhp5i^*£2 rafri Regenerates the Data and uses it to 

pJtfrnal % accuracv of ReconstI u-tion. 

Rnm1n^ I 2L T fi?°2? f I?t2tH!) , RESULT U 1) 


~ * * r i : ° A *' wAoyui , 1 ' w i / trtl I , HK Su L T I 

c ommon TH(in.sni) 

nSrllS !i8IJ*23' ACCESS* 'SEQl S' , FI 1.6* * TUFT*. ' ) 

d!EM4!'!H thci) ,i=i,in 

RtAD(23,* S(I),i=l,{l) 


r> ri % 'ri j a a J 7»f ^ v v *- ■ a c# v * 

thii) , 1 = 1 , 11 ) 
rlad( 23,*5( SCI), in, h) 


OUTPUT' ) 


ft=o,n ' v 
go io i=i,h 

B = COSDCTH(D) 

EPSPEL=0*01 
tPSABS*0,01 

? S sR d k >0 }r ^ i F / , A ' S » . ft ■ b ■ s » e 1 ps 1 r e i l , h , es : o ■ r, r a ) 
cnlliNUp' ,?00,XFAiI,) 

do' 29 1=1,11 
^SULTCnsZ.O^RESULl'Cn 

2 ?It^ c 5 JT 5 ??*^ 8 f If i >®f - 2 ? W TJ t« Pi **E-' l *»ci»AT- » 

STOP 

EfiD 

FUNCTION F { X) ' 

COMMON TH(tl) ,Stll) 

K=10*X 

iJl S f 1 ^ <j? \ 

F=S(J)+C10,*X-K)*(S(J+1)-S(J)) 

RETURN 



APPENDIX— E 


Referring to the figure we can write. 


J<Coi9 


3,k 


= 2 (x 2 +R 2 sin^ ) k dx 




= 2(x 2 +R 2 sin% ) k x I =k^ (x 2 +r 2 Sin%-) k " 1 (2x 2 )dx 

= 2 R 2k+1 Cos9 “2k /%^+ R 2 Sin 2 £ ) k dx 
+2kR 2 Sin 2 0 / (x 2 +R 2 Sin^ ) k “ 1 dx 


<4^ = 2R 2k+1 Cos “2k d q , k +2kR 2 Sin| dg" 


7 


k“l 


r 

= 1 2R 2k+1 Cosfc + 2kR 2 Sin%-d 

2k+1 l 


a,k“i 


= 2R Cos 

" 2k+l ” 


R 2 k”i + k sin Tan d 


,k“l 


For ease in FORTRAN programming. 


J k (i) = 2R 2 Cos0 j \ R 2k ~ 3 +(k-l)Sinfr ± Tan (k 

~2k“ l ' L J 


= 2R Cos 9 ± 


k= 2(1) - n 
i= 1(1) - n 
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